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ABSTRACT 

We present an analysis and interpretation of the Kepler binary system KOI 1224. This is the fourth binary 
found with Kepler that consists of a thermally bloated, hot white dwarf in a close orbit with a more or less 
normal star of spectral class A or F. As we show, KOI 1224 contains a white dwarf with Tgff = 14, 400 ± 1100 
K, mass = 0.20 ± O.O2M0, and radius = 0.103 ± 0.004 T?©, and an F-star companion of mass 1.59 ± 0.07 
that is somewhat beyond its terminal-age main sequence. The orbital period is quite short at 2.69802 days. The 
ingredients that are used in the analysis are the Kepler binary light curve, including the detection of the Doppler 
boosting effect; the NUV and FUV fluxes from the Galex images of this object; an estimate of the spectral type 
of the F-star companion; and evolutionary models of the companion designed to match its effective temperature 
and mean density. The light curve is modelled with a new code named I car us which we describe in detail. Its 
features include the full treatment of orbital phase-resolved spectroscopy, Doppler boosting, irradiation effects 
and transits/eclipses, which are particularly suited to irradiated eclipsing binaries. We interpret the KOI 1224 
system in terms of its likely evolutionary history. We infer that this type of system, containing a bloated hot 
white dwarf, is the direct descendant of an Algol-type binary. In spite of this basic understanding of the origin 
of KOI 1224, we discuss a number of problems associated with producing this type of system with this short 
of an short orbital period. 

Subject headings: binaries: eclipsing — stars: evolution — stars: individual (KOI 1224) — techniques: photo- 
metric — white dwarfs 



1. INTRODUCTION 

Close binary systems containing a single compact star (i.e., 
white dwarf, neutron star, or black hole) have the potential 
to enhance our knowledge of stellar evolution, perhaps even 
more so than for double degenerate systems whose evolution 
is yet more complicated. In almost all cases, the compact 
star originated from the remnant core of the original primary 
star in the binary system. The primary's envelope was lost 
via some combination of Roche-lobe overflow mass transfer, 
stellar wind, common envelope, or explosive process (e.g., a 
supernova). 

In cases where the original secondary star is now transfer- 
ring mass back to the compact star, the systems may be highly 
detectable via the release of gravitational potential energy in 
the form of optical, UV, or X-ray radiation. However, systems 
in which the compact star has not yet started to accrete may be 
much more difficult to discover. This is especially true in the 
case of white dwarfs in orbit with intri nsically much brighter 
stars that outshine them (e.g., Regulus;'Gi es et al.|2008t|Rap \ 
paport et al.|2009| ). Even in the case where the observer is m 
a fortuitous orientation to possibly see eclipses, such events 
may be of very small amplitude, e.g., at the 10~^ - 10~^ level 
if the white dwarf has cooled to near the temperature of the 
parent star. Due to the extraordinary photometric precision 
of the Kepler mission, however, eclipses and transits of this 
depth are readily detectable. 

Since the Kepler mission has been launched, there have in 
fact been a number of such binaries containing a white dwarf 
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that have been discovered (|Rowe et aT]|2010[ |van Kerkwijk| 
^raL]l2QT0l [Carter et al^|2011| ). InlEose systems the white 
dwart is still sufficiently hot (i.e., - 12000- 18000 K) and 
thermally bloated (i.e., up to ~ 10 times their degenerate 
radii) that the eclipses and transits are in the range of up to a 
couple of percent of the system light. We report on the study 
of a fourth such system with a nu mber of similarities to the 
others, KO I 1224 (KIC 6606653; [Borucki et al.||2QrT| [Prlal 
al.|2011|) ; however J t has the shortest among the orbital pe- 
)ds at 2.69802 dayqj and the primary star in this system has 
the lowest effective temperature and is the most evolved away 
from the main sequence of any of the previously detected hot 
white dwarf systems. 

We report here an analysis of the Kepler Ql and Q2 pub- 
lic data for KOI 1224, as well as supplementary observations 
from the GALEX satellite. We make use of stellar evolution 
models to better estimate the mass and evolutionary state of 
the current primary F star in the system. 

We analyze the Kepler light curve using Icarus, a newly 
developed state-of-the-art binary light curve synthesis code 
which is aimed at the study of detached and semi-detached 
binaries. This code has the advantage of handling spectro- 
scopic data as well as Doppler boosting, irradiation effects 
and transits/eclipses, hence making it particularly attractive 
for irradiated binaries. An eclipsing binary such as KOI 1224 
represents an ideal testbed to benchmark the synthesis code 
because its results can be easily compared with those of a 
semi-analytic analysis. The synthesis code is described in 
some detail in this work. 

The organization of the paper is as follows. In ^we intro- 
duce the new binary light curve synthesis code. Some detailed 

5 We note that a new binary system, ISWASP J024743.37-25 1549.2, re- 
cently discovered with the WASP survey (Maxted et al. 2011) contains a 
bloated hot white dwarf in an even shorter period binary with Porb =0.668 
days. 
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specific aspects of the code are given in Appendix |A| After 
presenting KOI 1224 and the available data in ^ we proceed 
with the light curve modelling in ^ first semi-analytically 
and then with our synthesis code. We discuss the proba- 
ble evolution of KOI 1224, its implications for understand- 
ing mass transfer, and its relation to Algol systems in ^ We 
summarize our findings and draw some general conclusions 
in^ 

2. BINARY LIGHT CURVE SYNTHESIS CODE 

We have implemented our own light curve synthesis code 
with the purpose of modelling the light curves and spectra 
of detached and semi-detached binaries. Our co de, called 
Icarua^ is mo st similar to the ELC code by Orosz & 
|HauschiIdt (2000). Essentially, the code constructs a tinite- 
element surface grid for a star having some pre-determined 
physical and binary properties by solving the gravitational 
equipotential equation. Each surface element is characterized 
by a set of physical properties (i.e., temperature and an effec- 
tive gravity), and the observed flux is obtained by integrating 
the specific intensity emerging from the surface visible by an 
observer located in a given direction. In the event that both 
stars contribute significantly to the total observed light, we 
model the second star of the system by inverting the mass ra- 
tio and adding a half orbital shift to it. 

The core implementation is made in Python and relies heav- 
ily on the Scipy and Numpy librarieq^ In order to bolster 
the execution speed, critical components are written in C and 
included directly in the Python code with the Scipy Weave 
modulq^ The modular, object-oriented approach that we 
have adopted should facilitate adding features to the code. 
Hence, the code currently includes various effects such as 
eclipses/transits, irradiation from a companion, ellipsoidal 
light variations, and Doppler boosting. As oppose d to other 
synthesis codes, such as the Wilson-Devinney code ( [Wilson & 



bevinney 1971), or its modern incarnation in PHOEBE (Prsa 



"Sl Z witter 2005), which bundle together the binary synthesis 



and the parameter optimization in a single piece of software, 
ours resembles more a suite of routines that can be glued to- 
gether in order to perform binary modelling. 

The main module of our code is the physics engine, which 
is responsible for synthesizing the star. A separate module 
deals with the atmosphere model and provides utility func- 
tions that return the intensity given a set of properties - typ- 
ically the temperature, surface gravity, emission angle and, 
sometimes, the velocity. This layer works independently of 
the binary light curve physics and serves only as a frontend 
to retrieve the specific intensities from a synthetic atmosphere 
model lookup table or an analytic model (e.g., blackbody). 
These atmosphere models can be integrated over a passband 
or used with their full spectral content. In the same spirit, the 
physics engine does not tamper with the returned values and 
leaves it to the user layer to alter them if necessary (e.g., to 
resample the spectrum or apply reddening). The user layer 
exists to connect the physics engine, the atmosphere models 
module, and the experimental data. It typically contains (1) a 
data reader function; (2) a utility function that takes user input 
parameters, recasts them into the physics engine input param- 
eter format, and post-processes the model data to match the 
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observed ones; and (3) a function that returns the fit residuals 
or goodness-of-fit. Finally, on top of the user layer, a fitting 
layer can be added. The function that returns the goodness- 
of-fit can be fed within the framework of a minimization pro- 
cedure, a brute-force search, a Bayesian optimization, or any 
other algorithm suited to the needs of the problem. 

The main innovation of our light curve synthesis code is that 
it can handle phase-resolved spectroscopic light curves and 
hence it allows one to perform modelling without the need for 
fitting radial velocities separately using a template and hav- 
ing to introduce a correction factor for the displacement of 
the light-center with respect to the barycenter. This capabil- 
ity will be demonstrated in a forthcoming paper. Thorough 
details about our code are provided in Appendix [A[ 

3. DATA AND DATA REDUCTION 

We have obtained the KOI 1224 data from the Kepler pub- 
lic archiv^ The original raw light curve, presented in Fig- 
ure fl] includes the first and second quarters (Ql and Q2) data, 
which were collected using a 30-min integration time. We re- 
moved the offset between the two quarters and corrected for 
three jumps in the flux. Next we cleaned the raw light curve 
from outliers using the following method. We removed the 
systematic trend using a 7th-degree polynomial and folded 
the light curve at the orbital period. We calculated the run- 
ning median and running standard deviation statistics using 
a 20-point window in order to identify outliers deviating by 
more than 1.5 standard deviations of the data scatter around 
the median, which we flagged and discarded for the rest of 
the analysis. This threshold was chosen in order to remove 
the obvious visible outliers as some data points were clearly 
beyond the regular data scatter. In total, the outliers removal 
discarded 699 out of the 5723 data points. 

3.1. Ultraviolet GALEX Data and Other Photometry 

In order to constrain our modelling, we started with infor - 
mation from the Kepler Input Catalog (IBrown et al. 201 1| ), 
which lists optical griz and infrared JHK 2MASS photometry 
as well as estimates of the temperature Tgff = 6352K and red- 
dening Eb-v = 0.120 (see Table [T]). The reddening is inferred 
using a simple exponential dust model, which, for this source, 
is consistent with the total dust column of Eb-v = 0.114 ex- 
pected along this line of sight from the COBE/DIRBE and 
IRAS based dust maps of'Schleg el et al.| ( [T998j ). Comparing 
the griz and JHK colours (dereddened using the coefficients 
of Schlegel et al. 1998 t to thos e inferred from spectrophotom- 
etry of MK standards ( [Covey et al._2007) , we infer a spectral 
type F6V, with an uncertainty of 1 subclass, corresponding to 
a temperature of 6500 ± 200 K ( |Cox|l20QQ| ). This is consis- 
tent with that of the Kepler Input Catalo g, within the expected 
uncertainty of ~200K (Brown et al. 2(311). From a spec- 
trum taken at the 1.6-m Observatoire du Mont-Megantic (L. 
Nelson, private communication) we infer a spectral type F5- 
6 III-IV (more likely luminosity class IV ), which implies a 
temperature similar to that quoted above ( [Gray et aL]|2001| ). 
Below, we will use Tiff = 6350 ± 200 K. 

We checked other archives for additional information, and 
found that KOI 1224 was detected by GALEX (see Table [T]), 
with ([FUV] - [NUV]) = 1 .41 ± 0.02 and ([NUV] -g) = 3.55± 
0.02. The estimated intrinsic ([NUV] - g)Q = 2.9 color is bluer 
than expected for the primary (^ 4.0; [Bianchi et al.||2007 
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Table 1 

Archival Photometry of KOI 1224 



Filter Magnitude^ 

"gaZeXFUV^ 18.98(2) 

GALEXmjy'' 17.57(1) 

KlCg^ 14.015(20) 

KlCr^ 13.651(20) 

KIC 13.589(20) 

KICz^ 13.556(20) 

2MASS7^ ... 12.77(2) 

2MASS//^ .. 12.56(2) 

2MASS ... 12.52(3) 



Color Reddening^ Intrinsic 



1.41 


-0.09 


1.50 


3.55 


0.61 


2.94 


0.364 


0.125 


0.239 


0.062 


0.080 


-0.018 


0.033 


0.073 


-0.040 


0.79 


0.15 


0.64 


0.21 


0.01 


0.20 


0.04 


0.01 


0.03 



The numbers in parentheses are the formal 1-cr confidence uncertainties 
and are applied to the last significant figures given. 

^Magnitudes are in the AB system for GALEX and KIC, and Vega for 
2MASS^ ^ 

I otherwise 



t Rev et al.H2007) for FUV, NUV,|Schlegel et al, 

^UALLA archive (http : //galex". srsci 



(http 
^Brown et al. (2o Tl| 
TSkrutskie et al. ( 2006 » 



1998 



■mi 



Vennes et al. 201 1 ), but we cannot exclude that it could signif- 
icantly contribute to the NUV flux given the uncertainties in 
the temperature, the reddening and the extinction curve. How- 
ever, no contribution from the primary is expected to the FUV 
flux, and indeed, with ([FUV] - [NUV])o = 1.50 ± 0.03[3the 
source is bluer in the UV than any of the brighter stellar 
sources near it. Thus, the FUV emission is almost certainly 
from the hot white dwarf secondary. 

In order to estimate a temperature for the white dwarf, we 
can use the fact that when the white dwarf is eclipsed, the 
flux in the Kepler bandpass decreases by 1.4%. Here, the 
1.4% drop includes a 10% correction for a "third light" (e.g., 
a fainter blended star within the Kepler postage stamp), and 
hence it represents the white dwarf's contribu tion to the sys - 
tem's luminosity. From color terms of Brown et al. ( 1201 1| ), 
the r-band flux should decrease about 10% less. We thus es- 
timate r 18.4 ± 0.2 and infer ([FUV] -r)c^ 0.6 ± 0.2 and 
([FUV] -ro)^0.0±0.3. We can compare this with predicted 



GALEX and optical colors for white dwarfs from | Vennes et al. 
(2011 ). For their lowest mass, O.4M0, they Hst ([FUV]-r) = 
3.29,1.90,0.86,0.27,-0.03,-0.26,-0.46 for Teff = 10000 to 
16000 K in steps of 1000 K, with the color depending some- 
what on gravity below 13000 K (here, we converted the V- 
band magnitudes of Ve nnes et al.| ( |2011| ) to r- band using V - 
r - 0. 19, valid for temperatures m this range; [Fukugita et al. 
1996]). Thus, we infer a temperature of - 14000 zb 1000 K. 

For the above temperature, the models also predict 
([FUV] - [NUV])o = 0.0 ± 0.1, which is much bluer than the 
observed color. Apparently, therefore, either our estimate is 
inaccurate and the true temperature is 10500 ± 500 K, or the 
primary contributes ~75% of the NUV flux, i.e., is about a 
magnitude brighter than expected (perhaps because of stellar 
activity associated with rapid rotation). From the ratio of the 
transit to eclipse depths (see |4 1| ), it seems clear the latter is 
the correct interpretation, andtEe white dwarf effective tem- 
perature is ~ 14000 ± 1000 K. This i s fu rther confirmed by 
our formal fits to the light curve (see |4.2|). 



4. MODELLING THE LIGHT CURVE OF KOI 1224 
4.1. Semi- Analytic Modelling 



We first performed a semi-analytical analysis of the 
KOI 1224 light curve in order to obtain rough estimates of the 
system parameters. We fitted the raw data (see Figure [l]) for 
three out-of-eclipse orbital harmonics as well as a 7th-degree 
polynomial and three "pulsations" (at 3.49, 1.75 and 1.17 d). 
The polynomial function aims to remove the long-term in- 
strumental variations, while the pulsations were detected as 
significant signals in a power spectrum. The 3.49-day pulsa- 
tion is strongest and might reflect the rotation period of the 
primary, as seemed to be the case for KOI 74, where a pos- 
sible pulsation period of ~0.6d matched the rotation period 
inferred from its rotational velocity. If so, it would imply the 
star is rotating sub- synchronously, perhaps because the star 
is expanding during the course of its evolution. We then re- 
moved the systematic trend from the data and folded them at 
the orbital period. 

We measured various fiducial points on the white dwarf 
eclipse profile using a detailed plot of the light curve (Fig- 
ure [2]) in order to estimate the radius of the primary star in 
units of the semi-major axis, i.e., Ri/a. We found the phase 
of the first contact, the start of the full eclipse, the end of the 
full eclipse, and the last contact to be ^ -0.0471, ^2 ^ 
-0.0369, 03 - 0.0365 and (j)4 - 0.0466, respectively. From 
these values, measured in orbital cycles from the primary's 
inferior conjunction, we can infer Ri/a = 7T((j)3 -(j)i) = 0.263 
and Ri/a = 7r(04 - (j)2) = 0.262 (implicitly assuming / = 90°, 
since we have no other way to infer it from simple geometry). 

The first three out-of-eclipse orbital harmonics capture 
information about the system properties that can be de- 
rived from ellipsoidal light variations, irradiation effects, and 
Doppler boosting (see the folded out-of-eclipse light curve in 
Figure [3]). From our fit we found the fractional amplitudes 
Al = (2.542 ± 0.058) x lO'^, A2 = (37.675 ± 0.071) x lO'^ 
and A3 = (3.085 ±0.074) x 10""^ along with the correspond- 
ing phases 0i = (-6.25 ±0.48) x 10"^, 02 = (2.6 ± 1.1) x lO'^ 
and 03 = (1.5 ± 1.0) X 10~^. The phases are in orbital cycles 
measured from the ascending node of the white dwarf for 0i 2 
and from the white dwarf eclipse for 03. Because the ingress 
and egress time scales are comparable to that of the Kepler in- 
tegration time, it would be futile to make any inference about 



R2/a without a proper numerical treatment, as in §4.2 

In the simple, semi-analytic analysis that follows, we cor- 
rect amplitudes Ai and A2 up by a factor of ~10% to ac- 
count for the third light in the system, and Ai up by an addi- 
tional ~6% to account for the much smaller expected Doppler 
boosting effect of the whit e dwarf which decr eases the total 
amplitude of Ai (see, e.g.. Carter et al. 2011[ equation 11). 



and ^2- 



We label these "corrected" amplitudes A 

The first out-of-eclipse orbital harm onic arises from 
Doppler boosting (see |A.7[ and van Ke rkwijk et al.|[20TQ 



[Loeb & Gaudi|2003| for more details) and has an amplitude: 



^icos0i =/db- 



(1) 



From the color information of KOI 1224 presented in |3.1[ we 
estimate /db = 3.55 at 6350 K and logg = 4. Using the above 
expression we find that the harmonic amplitude A[ implies 
a projected velocity amplitude Ki :^ 23 ± 1 km/s, where the 
cited uncertainty includes 3% each from the amplitude and 
/dr. From Ki we find the mass function: 



For a stan dard extinction curve, A(FUV,NUV)/£5-y = (8.16,8.90) 
jRey et al.|2007| . Given the small difference, the intrinsic color is secure. 



Ml sin^ / 
(Ml +M2)^ 



^0.00344 ±0.00040 



(2) 
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Figure 1. Raw light curve of KOI 1224 from Kepler (black dots) with a best-fit solution (solid line). The dashed line shows the 7th-degree polynomial and the 
three pulsation harmonics used to remove the systematic trend (shifted upward for visibility reasons). The lower panel displays the fit residuals. 
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Figure 2. Light curve of KOI 1224 (black dots) along with a best-fit solution (solid line) folded at the orbital period of the system. The systematic trend and 
the third light contribution have been removed from the data, and the flux range is adjusted to display the primary and the secondary eclipses in the right and left 
panels, respectively. The lower panels show the fit residuals. 
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Figure 3. Light curve of KOI 1224 (black dots) along with a best-fit solution (soHd line) folded at the orbital period of the system. The systematic trend and the 
third light contribution have been removed from the data, and the flux range is adjusted to trim the eclipses in order to enhance the visibility of the ellipsoidal and 
Doppler boosting components. The lower panel displays the fit residuals. 
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The second orbital harmonic is due to the ellipsoidal light 
variations, and has an expected amplitude: 



^2COS02=/eV 



M2 

m7 



sin^ /, 



(3) 



where /ev is a pre-factor of order unity that depends on the 
limb and gravity darkening: 



45 + 3w ^ 
20(3 -w) 



(4) 



( Morris||l9 85'p] Using the corrected amplitude, A'2 and the 
geometrically measured value of Ri/a, and considering that 
sin/'^ 1, we find a mass ratio ^ = M2 /Mi = (0.23±0.02)//ev. 
We also find that /ev — 1.59 if we take the linear limb dark- 
ling w ~ 0.56 for Kepler's bandpass at 6350 K from Sing 
(|2 010t, an d r ~ 0.66 by integrating Equation 10 from Mor- 
)| ( |1985| ) over the Kepler bandpass. Hence, we estimate that 
q = 0.144 zb 0.020. One also sees that the ellipsoidal term is 
phased as expected with the ascending node. 

Combining the mass function found from the A[ amplitude 
and the mass ratio from the A2 term, we find 



Ml 1.5 Mo M2::^O.22M0. 



(5) 



Second order terms from ellipsoidal light variations arise 
because the nose of the star (i.e., near the Li point) is darker 
than its back side. This results in the third harmonic, with a 
relative amplitude A3/A2 = 0.082 ±0.002. This values dif- 
fers slightly w ith that of 0.06 calculated with the equations of 
Morris! ( p"985| ), assuming the same coefficient as above. 

Lastly, there should also be a contribution to the fundamen- 
tal from the ellipsoidal light variation effect, with an ampli- 
tude of 3^^/5 = (2.04 ±0.05) x lO'^, with a maximum at 
the white dwarf's eclipse. Instead, the fundamental modu- 
lation is a bit shifted towards the white dwarf transit, with 
amplitude A[ sin((/)i) = (1.07 ± 0.08) x 10""^. This indicates 
that the shift of the fundamental modulation is mainly due 
to irradiation, rather than second-order ellipsoidal light varia- 
tions. Assuming the second-order contribution due to the el- 
lipsoidal light variation is correct, we infer an irradiation term 
of /irr (3. 1 1 ± 0.09) X 10~^. According to harmonic decom- 
position of the reflection effects in close binaries by Kopal] 
(fl959), we can write the fractional amplitude of the first har- 
monic term due to reflection as: 



fin 



•ai- 



L2 



11/^ 
3 4 V ^ 



(6) 



where a is the primary's albedo, L2/L is the ratio of the 
secondary's luminosity to the total luminosity. In this esti- 
mate we have ignored the leading-order {R2 / of term (due to 
the irradiation of the secondary by the primary) because it is 
nearly two orders of magnitude smaller than the L2/L(Ri/df' 
term - as we find in |4.2| (see Table [2]). From this we in- 
fer L2/L ^ 0.01 la, which is quite consistent with the ^1A% 
drop in the light curve during the white dwarf eclipse. 

In the next section we discuss our numerical modelling 
from which we obtain more accurate results for the system 
masses. The main reason for this improved accuracy is that 
we utilze measured properties of the primary star to directly 
determine its mass. In that case, the use of A2 in determining 

^ ^ Here we dropped the factor k from the original equation because it is 
negligible 



the mass ratio is less critical, and, more importantly, M2 then 
depends principally on the cube root of the mass function, and 
so is only linearly dependent on A[ (rather than on its cube). 

4.2. Numerical Modelling 

The numerical modelling has been accomplished using our 
new binary light curve synthesis code I car us sketched in ^ 
and described more completely in Appendix [A[ 

In order to fully describe the light curves, our code includes 
the following list of free parameters (with indices 1 and 2 re- 
ferring to the F-star primary and the white dwarf secondary, 
respectively), which are also summarized in Table|2] the mass 
ratio, q = M2/M1; the filling factoiF] /i,2 = Xnose(i,2)ALi(i.2); 
the (polar) temperature ratio, T2 / Ti - me primary's (polar) tem- 
perature, 71; the irradiation temperature ratio, Ti^i^/Ti (see 



§ A.4 for the definition of 7i irr); the orbital inclination, /; and 
a third light contribution, L3 . 



Table 2 

KOI 1224 Parameters 

Parameter Value 

Orbital period, Porb (d) 2.69802(2) 

Epoch of primary's inferior conjunction, Tq (MJD) 55030.038(1) 

Model Parameters 

Mass ratio, g = M2 /Ml 0.128(10) 

Filling factor, /i (fraction ofxn^) 0.389(6) 

Filling factor, /2 (fraction of xl2^) 0.034(2) 

Temperature, Ti (K) 6250(200) 

Temperature ratio, 72/ 7i 2.3(1) 

Temperature ratio. Tin-/ 71 0.164(3) 

Orbital inclination, i (degree) 85(1) 

Third-light, L3 (fraction of total light) 0. 10(5) 

Extra noise, (Xextra (fraction of average flux errors) 3.69(4) 

Derived Parameters 

Massprimary, Ml (M0) 1.59(7) 

Mass secondary, M2 (M0) 0.20(2) 

Semi-major axis, a (Rq) 9.9(2) 

Projected velocity amplitude, Ki (km s~^ ) 2 1 .0( 1 .7) 

Temperature, 72 (K) 14400(1100) 

Relative volume-averaged radius, {R\)/a 0.270(5) 

Relative volume-averaged radius, (/?2}/« 0.0103(4) 

Volume-averaged radius, {Ri) 2.67(6) 

Volume-averaged radius, {R2) 0.103(4) 

Luminosity, Li (Lq) 9.8(1.3) 

Luminosity, L2 (Lq) 0.41(13) 

Volume-average density, {pi> (gcm-3) 0.118(6) 

Volume-average density, {p2> (gcm-3) 270(30) 

Polar surface gravity, log gi (log J Q (cm s~^)) 3.79(2) 

Polar surface gravity, log gi (log 10 (cm s~^)) 5.72(5) 



The numbers in parentheses are the formal l-a confidence uncertainties 
and are appl ied to the last significant figures given. 
^See §A.2|for definition. 



One more parameter is required in order to fully specify the 
system dynamics and solve for the masses. This is required 
in order to determine the projected velocity amplitude of the 
primary, Ki , which, in turn, allows us to calculate the Doppler 
boosting amplitude. Instead of introducing this quantity as an 
additional free parameter to the fit, we use the primary mass. 
Ml, itself; this is possible because the mean stellar density 



See { A.2 for definition. 
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tential surfaces. We also used a gravity darkening coeffi- 
cient jg H^rv = 0.175 for both stars based on empi rical values 



from Che et al.| ^ll'p] As we pointed out in §4.1 irradi 



9000 8000 7000 6000 5000 

Effective Temperature (K) 

Figure 4. Evolution of the primary star in the 7i eff- (pi) plane. The 10 
curves are for stars starting on the ZAMS for masses between 1.1 and 2.0 
Mq in steps of 0. 1 Mq . The colored region represents a nearly continuous 
mapping of primary mass in the 7i eff - (pi) plane. This map was generated 
from the evolution tracks using an assumed 2-dimensional Gaussian probabil- 
ity distribution with l-a uncertainties of 200 K in 7i gff and 10% in (pi } . The 
black dot indicates the position of KOI 1224 from our numerical modelling 
along with the 68% confidence intervals on the temperature and density. 

of the primary, (pi), and its effective temperature Tefjo] de- 
termine a nearly unique mass and age, provided that the star 
has not evolved far up the giant branch. The average density 
of the primary star, in turn, can be found using Kepler's third 
law, which depends only on known quantities such as the vol- 
ume averaged stellar radius (7?i)/ap^ 



(pi) = 



(7) 



We have used the new ly developed MESA stellar evolution 
code ( [Paxton et al.||2011^ ) to evolve a sequence of stars be- 
tween lTr"an3"T^OM0 m steps of 0.1 Mq. These tracks are 
shown in FigureHin the (pi) -71 plane. Given that there is an 
uncertainty in both (pi) and 7i, and the fact that the evolution 
tracks cross each other near the TAMS, we have converted 
this discrete set of tracks to a nearly continuous mass distri- 
bution in the {pi)-Ti plane (see Figure [4]). To generate this 
finely gridded (400 x 600) mass array we generated a prob- 
ability weighting function associated with each point on the 
evolution tracks from a 2-dimensional Gaussian distribution 
in (pi) and 7i - as described in the caption to Figure Q If 
more than one probability was assigned to a given point in 
the (pi) - Ti plane for any given track, we took the highest of 
the probabilities. Finally, we used these probabilities to com- 
pute a weighted "mean" mass at each location in the (pi) - 7\ 
plane. 

We fixed the orbital ephemerides at the values reported 
in Table [2] and, because the rotational periods of the stars 
are unknown, we have assumed that they are rotating co- 
synchronously with the orbit for the calculation of the equipo- 

For this purpose, we approximated the effective temperature by 7i, the 
polar temperature. 

{Ri)/a is found after solving for the equipotential surface, which does 
not require knowing the masses in the system. 



ation effects on the white dwarf due to the primary are neg- 
ligible (~2% of the primary's irradiation by the secondary) 
and hence can be safely neglected. From exploratory fitting, 
we also realized that the flux uncertainties were insufficient 
to account for the root-mean- square of the fit residuals (i.e., 
the was large even though the fits visually looked good). 
Hence, we introduced an extra noise parameter, (Jextra. which 
was added in quadrature to the flux uncertainties. 

The process of calculating the goodness of fit of our 9- 
parameter model works as follows. We evaluate the syn- 
thetic light curve using a variable sampling scheme (i.e., more 
densely sampled near the transit/eclipse and less densely in 
the smoother ellipsoidal light variations region). This al- 
lows us to boost the computing speed while still capturing the 
critical details. The 30-min integration time of Kepler pro- 
duces a non-negligible smearing that needs to be accounted 
for. Hence, we resample the synthetic light curve at a higher 
resolution using a second-order spline in order to perform a 
30-min boxcar filtering. The synthetic light curve is then in- 
terpolated at the phase of each observed data point. Next, 
the systematic trend in the data is removed after fitting a 7th- 
order polynomial and three sinusoidal pulsations having un- 
known phases, amplitudes and harmonically related periods 
^puis., ^puis./2 and /puis./3, where the fundamental is ~3.49 d. 
Finally, the residuals are calculated using the detrended light 
curve. 

4.2.1. MCMC Algorithm 

We performed the parameter optimization using a Markov 
Chain Monte Carlo (MCMC) algorithm since it is well suited 
to making parameter inferences in high-dimensional prob- 
lems. Parameter correlations are sometimes complicated and 
hence we used a Gibbs sampler in our MCMC because it al- 
lows for an easier tuning of the pro posal distribu tion's accep- 
tance rate (see e.g. |Walsh 2004; G regory |[2005[ for more in- 
formation about MCMC and Gibbs sampling). 

We sampled the orbital inclination using flat priors in cos /, 
such that it eliminates projection effects. Based on t he Kepler 
Input Catalog estimation of the third light (see |3.1| ), which is 
listed as 10%, we used a logistic prior for this parameter: 

[l+exp((L3-0.15)/0.0l)]"^ . (8) 

In such way, the probability of L3 values in the range 0-0.1 
is close to unity before dropping to 50% at 0.15, and then be- 
comes negligible beyond 0.2. For the extra noise parameter, 
we used modified Jeffrey's priors [(Jextra+O - 1]. which provides 
equal weighting per decade. This choice is motivated by the 
scaling nature of the extra noise parameter. The constant of 
0.1 linearizes the prior when dextra becomes smaller than 10% 
of the average flux measurement errors. In a Bayesian param- 
eter estimation framework, the extra noise parameter has an 
effect similar to normalizing the reduced to unity in stan- 
dard fitting. Since Kepler uses a single broadband filter and 
the data contain no color information, we also made use of the 
photometry constraints (see |3.1| ) to add a Gaussian prior of 
6350 ± 200 K on the primary's temperature. Flat priors were 
assumed for the remaining parameters. 

This is appropriate for the primary star of 7i = 6350K, and the hotter 
white dwarf is sufficiently small that its gravity darkening has a negligible 
effect. 
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We ran a total of 6 independent MCMC chains, each con- 
sisting of 400,000 steps. A bum-in period of 20,000 steps 
was removed from each chain and we kept only every 360 
iterations - a process called thinning - in order to reduce 
auto-correlation effects. We compared the summary statis- 
tics of each MCMC to the others in order to ensure that 
convergence has been reached. For all model parameters, 
the G elman-Rubins convergence diagnostic ([Gelman & Rubin 



1992 ) yielded v ^ < 1.1, which indicates that the chains have 
reached stationary distributions. The final results presented in 
Table [2] and Figure |5| were obtained by merging all the chains 
together (after burn-m and thinning), while Figures |2] and [3] 
display an example of the best-fit light curve. 

We find constituent masses of 1.59 ± 0.07 M© and 0.20 ± 
0.02 M0 for the primary and white dwarf, respectively. The 
radius of the primary is 2.67 ± 0.067?©, indicating that it has 
evolved substantially off the zero-age main sequence, consis- 
tent with its low mean density of 0.12 g cm~^. The white 
dwarf's radius, 0.103 ±0.004 T?©, is ~5 times its degenerate 
radius, and indicates that it is still actively cooling (possibly 
punctuated with episodic shell flashes). In turn, this indicates 
that its cooling age is of order 120 - 600 Myr for masses in 
the range of 0.2 0- 0.22Mq (see, e.g., |Driebe et al.|1999[[Ner 



son et al.|2004 |. We discuss below what these properties may 
imply for the formation of KOI 1224. 

5. THE ORIGIN OF KOI 1224 

We have presented results for a system found with the 
Kepler mission that contains a thermally bloated hot white 
dwarf: KOI 1224 (KIC 6606653). The system consists of a 
normal F6 star that has evolved to near the TAMS in a 2.69 8- 
day binary with a white dwarf of mass O.2M0 and radius 
O.17?0. We discuss in this section how such a binary may 
have evolved. 

There are two plausible evolutionary paths to the forma- 
tion of the current KOI 1224 binary system. The transfer of 
mass from the white dwarf progenitor to the primordial sec- 
ondary was either dynamically unstable, leading to a common 
envelope phase, or it was stable, even if not completely con- 
servative. The former case leads to a dramatic shrinkage of 
the orbital separation while the latter results in only modest 
changes in the orbit. We consider each of these, in turn, and 
in more detail. For more discussion of the formation and evo- 
lution ofthese systems see Farmer & Agol ( 2003 ), Di St ^ano] 



( [2QTT] ), [van Kerkwijk et al.,(2010j and Carter et al., ( .201] 



5.1. Common Envelope Scenario 

In the common-envelope scenario, the primary star devel- 
ops a degenerate He core of mass '^O.2M0, starts to transfer 
mass to the secondary, and initiates a dynamically unstable 
process leading to a common-envelope phase. The secondary 
spirals into the common envelope of the primary, ejects the 
envelope, and this results in a very close orbit of the secondary 
with the He core of the primary (i.e., the current white dwarf 
of the system). The final orbital separation can be related to 
the initial (pre-CE) separation by 



McMs ri 



(9) 



(see, e.g. jde Kool|T99Ql|Podsiadlowski et al.l2003| ) where M^, 
Me, and Mp are the core, envelope, and total mass of the pri- 
mordial primary star, respectively, Ms is the mass of the sec- 
ondary, ri is the Roche-lobe radius of the primary in units of 



the orbital separation, and the product AacE encapsulates the 
energy ejection efficiency and binding energy of the common 
envelope. The latter CE parameter is often taken to b e ^1, but 
is more likely to < 0.1 for an evolved star (see, e.g.pewi & 



|Tauris|2000t|Tauris & Dewi|2001t |Podsiadlowski et a l.|2003| ). 
In this latter case, the orbit must shrink by a factor of > 100. 
For a final orbital period of 2.7 days, the initial period would 
have to have been > 5 years. Such a wide initial orbit for the 
progenitor binary would necessarily imply the development 
of a much more massive He (or likely CO) white dwarf of 
> 0.5 Mq. Thus, given the fact that the current white dwarf 
in KOI 1224 is much lower in mass, this seems to essentially 
rule out a common-envelope scenario. 

5.2. Stable Mass-Transfer Scenario 

Stable, but not necessarily conservative, mass transfer from 
the progenitor of the white dwarf to the secondary seems 
the much more likely route to the production of KOI 1224, 
though, as we discuss, this formation scenario is not without 
its own difficulties. In this scenario, when the primary fills its 
Roche lobe and commences mass transfer to the secondary, 
that mass transfer will proceed on the thermal timescale of 



the en velope of the do nor star (see, e.g., Podsiadlowski et al. 
|2QQ2i |Lin et al.|[2QTT] ) until well after the masses of the two 
stars have become equal. Thermal timescale mass transfer re- 
sults since the donor star is more massive than the accretor 
and it is not too evolved. The remainder of the mass transfer 
is driven by the nuclear evolution of the primary and is gen- 
erally m uch slower than the thermal timescale transfer (Lin] 
|etal.|20 11). 

Depending on the relative masses of the primary and sec- 
ondary at the time mass transfer commences, the outer en- 
velope of the accreting star may swell up to the point where 
further accretion is inhibited and a substantial fraction o f the 
transferred mass is ejected from the system (see, e.g. Kippen 



hahn & Meyer-Hofmeister|1977[|van Rensbergen et al.|2010 



201 1| ). If we define the mass retention fraction by the accretor 
to be /3, the value of (3 is likely to be low at the start of the 
mass transfer process and higher toward the end. 

As we have shown, the age of the white dwarf in KOI 1224 
is likely to be < 600 Myr bec ause it is still quite hot (see, e.g. 
Driebe et al. 1999; Nelson et al.|2004| ). The fact that the cur- 
rent primary star in KOI 1224 has a mass of ^1.6 Mq and 
has evolved to at least the TAMS (and probably somewhat 
beyond) implies an age of > 1.75 ± 0.3 Gyr for a single, iso- 
lated star of that mass. These two facts taken together imply 
that the current primary must have been already substantially 
evolved at the end of the mass-transfer phase when the white 
dwarf was unveiled. From all of this we can infer that the 
primordial secondary (now the current primary) was not sub- 
stantially less massive than the primordial primary. 

Therefore, while there likely was a brief phase of thermal 
timescale mass transfer, the bulk of the mass transfer should 
have been at a slower rate and therefore could have been 
mostly conservative. For purposes of simplicity and for pre- 
senting illustrative examples, we therefore parameterize the 
evolution of KOI 1224 as having a constant mass retention 
fraction, (3. We further take the specific angular momentum 
of any ejected material to be a in units of that of the binary 
system. 

Figure [6] shows how the period of such a binary evolves as 
mass is transferred from the primary, for a range of plausible 
values of /3. In all cases, the orbit shrinks as mass transfer 
gets under way. However, for /3 > 0.7 the orbital period even- 
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Figure 5. Results of the MCMC parameter estimation for the KOI 1224 light curve. The data from 6 independent chains were combined after applying a burn-in 
of 20,000 steps and a thinning of 360 iterations. Note that K\, Ri/a, R'lja, M\ and M2 are all inferred values from other model parameters (see ^4]for more 
details). 
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tually dramatically increases above the original orbital period. 
Whereas, for smaller values of /3 (< 0.5) the orbit either does 
not grow or continues to shrink with mass transfer. Regardless 
of the evolutionary path, the orbit cannot get close enough so 
that the accreting secondary overflows its Roche-lobe, other- 
wise the two stars would be in danger of merging. 
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Figure 6. Evolution of the orbital period of KOI 1224 for different con- 
stant values of the mass retention fraction, (3. In these evolutionary tracks, 
we assumed a single value of the angular momentum loss parameter, a = l. 
For illustrative purposes we chose the primordial masses to be \3Mq and 
1 . 1 Mq . For /3 ~ 0.5, the orbital period remains constant to within a factor of 
two over the entire evolution. However, for ^ > 0.6, the orbital period dra- 
matically increases during the later portions of the binary evolution. Under 
the assumption that a and /5 remain constant during the evolution, the mass 
transfer has to be relatively non-conservative, with /3 probably no larger than 
~ 0.6, in order to account for the short observed 2.7-day orbital period of 
KOI 1224. 



The range of possible progenitor masses is explored in Fig- 
ure [7] The shaded region indicates the most plausible pairs of 
progenitor masses, Mp and M^. The constraints leading to this 
region are that (i) the ratio of nuclear and thermal timescales 
of the primordial pair lies between 0.5 and 0.9; (ii) the ini- 
tial orbital period was shorter than ~10 days; and (iii) the 
primary expanded by at least 30% in radius by the time that 
Roche lobe overflow commenced. The red curves are con- 
tours of constant orbital period at the start of mass transfer, 
while the blue curves are contours of constant mass retention 
fraction, (3. We again assumed a single value for the angular 
momentum loss parameter of a = 1 . 

The rationale for constraint (i) above is that the secondary 
should have a nuclear evolution time that is not too much 
longer than that of the primary in order for it to have evolved 
off the main sequence only after the mass transfer was com- 
pleted, but before the white dwarf had a chance to cool below 
~ 14,000 K. The requirement (ii) for an initially short orbital 
period (^10 days) is to avoid developing a white dwarf mass 
that is too large (i.e., remaining below ~ 0.23 M©, our upper 
limit on M^d) to allow its progenitor to fit within its orbit at the 
time mass transfer commenced (see the discussion below). Fi- 
nally, constraint (iii) allows room in the initial binary for the 
primary star to develop a He core before thermal timescale 
mass transfer reduces the primary's mass below ^IMq and 
impedes its further evolution. 

Finally, in regard to the evolutionary history of KOI 1224, 
we note that there is a theoretical relationship between the 
mass of the white dwarf and the orbital period in systems 
where the progenitor of th e white dwarf was a lo^y-mass star 
(i.e., < 2.2Mq\ see, e.g., "Rappapor t et aP 1995 Lin et al. 
|2011| ) and where the mass transfer was stable. This is conve- 



^"1.0 




Figure 7. Diagram showing the initial progenitor masses of a binary system 
that could evolve to become KOI 1224. Plotted are the initial primary (pro- 
genitor of the WD) mass, Mp, and the initial secondary (now the ~ \.6Mq 
primary) mass, Ms. The red diagonal straight lines labeled 1 -0.3 (with a 
large font), are contours of the mass retention fraction, (3 (i.e., the fraction 
of the mass leaving the primary that goes to the secondary and is not ejected 
from the system). The curved blue contours with values ranging from 0.5 
to 10 (smaller font) are contours of constant initial orbital period (in days). 
The central shaded region provides the boundaries in the Mp -Ms plane that 
would yield the observed masses set by constraints on the ratio of thermal and 
nuclear timescales and the requirement that the primordial primary underfills 
its initial Roche lobe by at least a factor of 1.3 (to allow it to evolve before 
transferring mass). The evolutionary tracks were calculated assuming that 
the specific angular momentum of the ejected mass in units of the angular 
momentum of the binary per unit of the binary reduced mass is o; = 1 . 



niently summarized in Figure 5 of |Lin et al. ( 2011| ) and their 
equation (1). In their Figure 5 we can see that stars with ini- 
tial mass < 2.2M0 (the red and green points in that plot) pro- 
duce a steep functional relation between orbital period and 
the final core mass in the evolution of binary systems involv- 
ing neutron star accretors. This same relation also somewhat 
holds for donor stars of initial mass up to ~4 Mq that com- 
mence mass transfer in so-called "late case A" or "case AB" 
mass transfer, i.e., when the donor is near the end of its main- 
sequence phase. The fact that the accretors in the current 
problem are main sequence stars does not alter the core mass- 
radius relation on which this effect is based. 

If we utilize the expression relating orbital period to core 
mass: 

Porh ^ '—jj^ days (10) 



(1+25M3-5-F29M6) 



I Lin et al.|2011| ) we can infer that a final orbital period of 2.7 
ays should correspond to a white dwarf mass of O.2O6M0 
with a spread in the theoretical values of ~ ±0.015 M©. This 
is highly consistent with the mass of the white dwarf we find 
for KOI 1224. One caveat, however, is that the Porb-^wd 
relation i s based on the core mass-radius relati on for giants 
(see, e.g. Han et al.|1994t Rappaport et al.|1995] ). In order for 
the latter to hold, the mass transfer must be effected in at least 
a semi-controlled way (see 

6. CONCLUSION 

In this paper we presented KOI 1224 (KIC 6606653), a new 
Kepler eclipsing system consisting of a bloated white dwarf 
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having a normal, slightly evolved companion star in a com- 
pact 2.698-d orbit. Archival analysis of the optical and near- 
infrared photometry, as well as inspection of a preliminary 
classification spectrum (L. Nelson, private communication), 
reveal that the primary is likely an F6 spectral class ~ 6350^ 
star, while the GALEX UV flux allows us to infer a white 
dwarf temperature of ~ 14,000^, a value that is compara- 
ble with that found by our numerical modelling when using 
the primary's temperature as a prior. 

We modelled the light curve of this binary with Icarus, 
a new light curve synthesis code, which uses proper atmo- 
sphere models and a physical description of the binary pa- 
rameters that allows us to reproduce various effects visible in 
the light curve such as irradiation, ellipsoidal light variations, 
and Doppler boosting. The derived binary parameters agree 
relatively well with those inferred from a harmonic analysis 
of the Kepler data using semi-analytic approximations of the 
above effects. 

From our light curve fitting, we found masses of 1 .59 ± 0.07 
and 0.20 ± O.O2M0 for the F-star primary and the white dwarf 
secondary, respectively, in a 2. 69 8 -day orbit. These values 
pose a modest challenge in terms of binary evolution. Indeed, 
the low mass of the white dwarf likely excludes a common 
envelope scenario, which would require a more evolved pro- 
genitor - and therefore a more massive white dwarf - in or- 
der to produce a compact system with the observed param- 
eters. At the same time, the orbital separation of a system 
like KOI 1224 is expected to increase significantly in a stable 
Roche-lobe overflow evolution unless the fraction of mass ac- 
creted by the current primary remains relatively low, /3 < 0.5. 
In any case, it appears that the system had to start its mass 
transfer evolution at an orbital period that is not longer than 
~10 days (as discussed above). 

KOI 1224 is the 4th such system found with Kepler. The 
orbital periods are: 2.7 d (KOI 1224), 3.3 d (K HWD3), 5.2 
d (KOI 74), and 24 d (KOI 81) (current work; [Carter et al.| 



2011[ I van Kerkwijk et aL]|2010| ). The corresponding white 
dwarf masses are ^0.21, 0.26, 0.21, and 0.3 M©, respectively, 
with 10-20% uncertainties. Unfortunately, the theoretical re- 
lationship between orbital period and white dwarf mass (dis- 
cussed above) is too steep, i.e., Porb ^ ^wd' where ^ > 7, to 
test effectively on such a relatively narrow range of Porb and 

Preferably, one would like to find more such systems in 
wider orbital periods in the Kepler data where (i) the white 
dwarf masses are expected to deviate more substantially from 
~0.2 M0, and (ii) the Porb-^wd relation should work more 
robustly because of the higher likelihoo d of purely stab l e mass 
transfer. In this regard, we note that |Maxted et al.| ( |2011| ) 
have recently discovered a possibly related system with an or- 
bital period of only 0.688 days, and a highly bloated white 
dwarf of 0.29 ± 0.005 M©, as well as several other similar 
very short period systems (P. Maxted, private communica- 
tion). The Porb -^wd relation is difficult to quantify theoreti- 
cally at these very short periods (of < 1 day) because it is not 
yet clear whether these systems result from stable mass trans- 
fer, a common envelope scenario, or some less catastrophic 
form of unstable mass transfer. The Porb -^wd relation works 
only when mass is removed from the progenitor of the white 
dwarf in a reasonably controlled way such that the core-mass- 
radius relation, on which it is based, remains valid. 
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APPENDIX 

ICARUS: BINARY LIGHT CURVE SYNTHESIS CODE 
Stellar Surface Grid 

The underlying stellar surface grid is constructed u sing the triangular tessellation of a unit sphere that was inspired by the 
GDDSYN synthesis code ( [Hendry & Mochnacki|1992 |). The base vertices are initially generated from the primitives of an icosa- 
hedron. Each triangular surface is then sub-divided into four triangles by creating a new vertex at the midpoint of each edge. 
The new vertices are finally projected back onto the unit sphere by re-normalizing them to unity. This sub-division process can 
be repeated iteratively in order to obtain the desired number of surface elements. At a sub-division level n (where ^ = is the 
base icosahedron, and ^ = 1 is one sub-division further), the number of surface elements is ^faces = 20 (4^^), while the number of 
independent vertices is ^vertices = 2+10(4"). The coordinates of each face is calculated as the centroid of the triangle. 

There are several advantages to using the algorithm above. First, a triangular tessellation does not suffer from the uneven surface 
sampling problem common to simple meshes derived from an equally spaced grid in spherical coordinates. To circumvent this 
problem, other synthesis codes, like PHOEBE, use a variable sampling strategy with fewer points near the poles. The subdivision 
algorithm of a triangular tessellation is relatively easy to implement and, using an icosahedron as the primitive, yields triangular 
faces that are organized on a regular tile of hexagons and pentagons, and that have fairly equal areas. Triangular faces are by 
definition convex and hence are simpler to work with when it comes to dealing with transits and eclipses. For reasons that will 
become obvious later we also keep track of the associativities: each surface is associated to three vertices, and each vertex is 
associated to either 5 or 6 faces. 



Equipotential Surface 

As in the ELC code, our stellar surface is defined by the gravitational equipotential equation from |Avni & Bahcall| ( |1975| ), 
which takes into account the effects of stellar rotation: 



GM, 



1.1 

n r2 



-qx-\- 



{r 



(Al) 



where ri and r2 is the distance of a point measured from the barycenter of the star that is being modelled (labeled 1) and from 
the barycenter of its companion (labeled 2) in units of orbital separation, ^ = M2/M1 is the mass ratio, and ft = c^star/^orbit is the 
co-rotation factor expressed as the ratio of the stellar to orbital frequency. We work in a coordinate system having its origin at 
the barycenter of star 1, in which the x-axis points towards star 2's barycenter, the z-axis is along the orbital angular momentum, 
and the j-axis is along the orbital plane orthogonal to the x and z axes. Our code is currently limited to circular orbits but could 
easily be extended to elliptical orbits as in PHOEBE, for example. 

Initially, the position of the LI point, xli, along the x-axis is f ound by identifying the saddle-point of the above equation 
between stars 1 and 2. Then, as explained in |Qrosz & Hauschildt| ([2Q00 ), the potential corresponding to the "nose" of the star 
(the inner point towards star 2) is defined by specifying a tilling factor / = XnoseALi, which is expressed as a fraction of the LI 
distance. From there, the radius of star 1 can be evaluated in the direction of every vertex and face of the tessellated surface. 

The effective surface gravity is also calculated for each face using the equation g = || VV^|| as well as the components of the 
normal to the local surface. 



Surface Temperature and Gravity Darkening 
We assign a fiducial temperature to each face that takes into account the gravity darkening according to the equation: 



T{x,y,z) 

^ole 



(^pole 



(A2) 



where /3 typically varies between 08 and 0.25 depending on w hether the star ha s a convective or a radiative envelope, respec- 
tively ( |Lucy 1 1 967} | von Zeipelj 1 924] for empirical constraints, see |Che et al.|2011| ). 



Irradiation from a Companion 

We treat the effect of irradiation from a companion as a source of energy that increases the temperature of the star while 
preserving thermal equilibrium (i.e., no thermal inversion that gives rises to emission lines). In this case, irradiation changes the 
temperature distribution of the star as following: 



T'{x,y,z) = 



T{x,y,zf + n, 



cosx 



1/4 



(A3) 



f >eech||1985| ) where cosx is the angle between the surface normal and the direction to the companion (i.e., star 2), and Ti^ 
e irradiation temperature which describes the temperature increase. Note that working with T^r is more convenient than usir 



- IS 

temperature which describes the temperature increase. Note that working with T^^ is more convenient than using 
an irradiation luminosity and a bolometric albedo since empirically one effectively measures the back and front temperatures 
of a star using multi-band light curves. One can convert the irradiation temperature into an irradiation luminosity after the fact 
using Lirr = AiTcra^T^, where cr is the Stefan-Boltzmann constant and a the orbital separation, and work out the stellar albedo if 
the companion's luminosity is known such that a = Lirr/L2. If one prefers, it is also possible to fix an albedo and perform the 
calculation the other way around in order to provide an irradiation temperature to the code. 
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Our treatment of irradiation is different than that of PHOEBE and GDDSYN, which calculate the reflection effect of each surface 
element due to the other star in an iterative way. While their technique is certainly more accurate, it comes at the cost of a 
substantially larger computational burden and only becomes more relevant in the case of contact binaries, where one needs to 
account for shadowing effects due to the bridge of matter in order to be self-consistent. Moreover, the fact that heat might be 
partly redistributed over the surface in a way that is non-trivial to calculate from first principles helps justify our use of this 
approximation. 



Surface Area 

We pre-calculate the area of each surface triangle on the unit sphere. Once the equipotential equation has been solved we 
rescale each pre-calculated value by its proper radius in order to obtain the effective area of the surface elements. By this means, 
the integrated flux simply corresponds to a discrete summation over the visible surface. 



Atmosphere Models 

Once the stellar grid has been constructed, and the effective temperature and surface gravity determined for all the faces, we 
can evaluate the flux perceived by an observer located at a given orbital inclination and orbital phase. The backend of the code 
that returns the flux for each face can be swapped between a blackbody or an atmosphere grid. 

For the purpose of this paper, we have worked with the BTSettl atmosphere models of [AUard et al.| ( [2003 [ [2007 [ |2010 1) which 
are available online on the Phoenix web simulatoip^ We used these atmosphere models to cover a grid over the range 1000- 
15000 K in temperature and 3.0-4.5 in \ogg. Since only integrated spectra were available, we used the empirical limb darkening 
relationship of Neckel (2005 ). 

Prior to performing the modelling of the KOI 1224 data, we integrated the spectral grid over the Kepler passband. Such 
integration can be performed for any photometric filter that one wishes to use to model the light curves. However, our code 
can also handle working with a full set of spectral data. In which case it is possible to model full orbital-resolved spectroscopic 
data. When this is the case, the spectrum of each face is Doppler shifted at the appropriate velocity and hence the rotational 
broadening, orbital Doppler shift, and displacement of the light-center away from the bary center due to irradiation naturally arise 
from our model spectra. In princ iple, the R ossiter-McLaughlin effect (Rossiter 1924; McLaughlin 1924 ) (as well as its Doppler 
boosting analogue ( van Kerkwijk et al.|2010^,Shporer et al.|201 1) ) would also appear from spectra modelled in combination with 
the transit calculation described below. 



Doppler Boosting 

Doppler boosting is another feature implemented in our code. For most cases, the amplitude of the effect is negligible, and 
hence it is turned off to speed up the light curve calculation. For the KOI 1224 data, however, Doppler boosting is clearly detected 
and hence it was i ncluded in the following w ay. 

As described in | van Kerkwijk et aLlp010| ), Doppler boosting can be written as a modulation of the observed flux: 

/=l-/DBsin(^v/c, (A4) 

where (j) is the orbital phase measured at the inferior conjunction of the star, v is the orbital velocity, c is the speed of light, and 
/db is a factor that depends on the spectrum o f the source and the observing passband. Based on our analytical estimate of the 
primary's temperature and surface gravity (see |3.1| ), we evaluated the coefficient /db at temperatures of 6200 and 6500 K for \ogg 
of 3.5 and 4.0. The coefficient is a rather smooth function of the temperature of surface gravity, and hence we simply performed 
a bilinear interpolation to fetch the coefficient to apply to each face of the stellar grid before the total flux was integrated. In the 
future, we are planning to include the exact calculation at each temperature and \ogg value. 



Eclipse and Transit Calculation 

Eclipse and transit calculation are also included in our synthesis code. Each time that the stellar surface grid is recomputed 
with new input parameters, we perform several tests to evaluate whether partial or total occultations are possible as well as the 
orbital phase range within which they can occur. This allows for our code to revert to algorithms optimized for each situation. 

At an orbital phase (j) and an inclination /, the projected distance between the bary center of the stars is given by {5 of' = 
a^ (cos^ /cos^ + sin^ /) , where ^ = at the inferior conjunction of the primary. 

First, we test whether any kind of occultation can occur. We consider the maximum dimension of each star in order to obtain 
a conservative estimate. Hence, for any orbital range where Sa < ri,max + ^2,max, one star might be occultated by the other. If this 
is the case, a second test aims to check if a full occultation occurs, and, if so, what orbital range does it cover. If ri max < ^2,min, 
star 1 will definitely get fully eclipsed in the orbital range that meets the requirement (5a + ri max < ^2,min. The same test is also 
performed with the indices reversed. 

With these criteria identified, the flux calculation can be classified into three categories: (1) no occultation, (2) partial occul- 
tation, and (3) full occultation. Case 3 is the simplest as the star does not contribute to the total observed flux and hence no 
calculation is required. Case 1 is also relatively simple as every face that is not beyond the stellar limb will contribute; a standard 
flux integration over the visible surface of each star is performed. Case 2 is more complex and we explain our algorithm below. 



Available at 



http : //phoenix . ens-lyon . f r/ simulator/ 
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Figure 8. Diagram illustrating an occultation of the primary (larger star) by the secondary (smaller, blue star). The primary's tessellated surface is displayed with 
the following details: the triangle symbols mark the center of the triangular surfaces, the black dots are the vertices of the triangular surfaces and the red dots are 
the vertices of nested triangular surface sub-divisions. To evaluate the visibility of the surface elements, the outline of the secondary is found and a lookup table 
r^iO) created, where Q and are the polar coordinates of the limb in the sky plane, using the secondary's barycenter as origin. The projected distance of each 
vertex from the secondary's barycenter is then calculated. Surface elements having all their vertices lying within the light blue region are not visible. For those 
near the boundary, a weighting is applying. For the surface element highlighted in grey, one of its primary vertices is occulted and so is one of its sub-vertices. 
Hence, 2/6 vertices are hidden which implies that this surface element only contributes 33% its normal flux. 

Partial Visibility 

Evaluating the partial visibility of a stellar surface occultated by an other is a computationally intensive task that can be tackled 
in several ways. Our algorithm is meant to provide a balance between accuracy, execution speed, and ease of implementation. 

In short summary, our algorithm aims to calculate the fraction of a surface element that is visible to an observer. The first step 
consists in the determination of the outline of the star located in front (here we choose star 2 for the explanation). We make the 
approximation that the orbital range within which the occultation of star 1 occurs is small enough that the outline at the onset of 
the occultation is the same as at conjunction. This simplifies the calculation considerably since, when the orbital inclination is 
sufficiently close to edge-oif^ the limb of star 2 corresponds to its outline in the y-z plane at x = 0, which is easy to compute. 
In this fashion, we obtain a lookup table of r2(0) values, where and r2 are the polar coordinates of the limb of star 2 in the sky 
plane calculated using its barycenter as the origin. 

In the second phase, we compute the projected position of star I's vertices, 6i and ri, with respect to star 2's barycenter. If a 
vertex of star 1 lies within the limb of star 2 - that is ri < r2(0i) - then it is hidden. To calculate the partial visibility of a surface 
element, we make the approximation the visible fraction of the surface is equal to the fraction of its vertices that are not hidden. 
Since each individual vertex is associated with either five or six faces, we can use the vertex-to-face associativity to assign the 
partial visibility weighting factor to all the faces without having to perform redundant calculations. 

Our partial visibility algorithm relies on a simple projected distance determination for all the vertices of the visible face. When 
such a vertex is occulted, an extra weighting factor is added to five or six array elements corresponding to the associated surface 
elements. The computational cost of this technique remains small compared to a more accurate calculation such as that performed 
by GDDSYN. Moreover it can easily be parallelized. 

The main caveat of our algorithm lies in the limited spatial resolution of the grid faces; the partial visibility goes in 1/3 
increments. Hence, the light curve can display important jitter if the chosen grid resolution is too coarse. A notable symptomatic 
case is that of stars having very unequal relative sizes and for which the projected area of the smaller star becomes comparable to 
the area of a surface elements of the larger star. The obvious solution relies on using a higher resolution for the larger star, though 
this comes at the expense of an increased computational burden. 

To work around these problems, we have designed a method that takes advantage of the sub-division tessellation algorithm. 
Since the sub-division to a finer level adds vertices at the mid-points of the surface sides, it is easy to keep track of the associations 
of these vertices with the parent surface elements at the lower resolution. Using such a nested sub-division yields partial visibili- 
ties that can be computed to a much better precision than the original 1/3 step-size, increasing to 1 /6 for one extra sub-division 
and 1/15 for two. This algorithm comes with a relatively small computing footprint: since the surface properties (temperature 
and gravity) vary smoothly, it is sufficient to obtain an accurate partial visibility weighting for a low-resolution grid using the 
sub-division. In other words, there is little if no gain in calculating the emergent intensity for all the surface elements at the higher 
resolution. 



17 A more accurate outline could be calculated but our approximation is ^^^^ ^^^^^^ ^^^^ ^224 



